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ABSTRACT 

This report shows, in detail, how the geomagnetic field interacts with the particle flux of the radia- 
tion belts to create a hazard to spacecraft and humans in near-Earth orbit. It illustrates the geometry 
of the geomagnetic field lines, especially around the area where the field strength is anomalously 
low in the South Atlantic Ocean. It discusses how the field will probably change in the future and the 
consequences that may have on hazards in near space. 


iii 


The Geomagnetic Field... 




CONTENTS 


LIST OF FIGURES vii 

LIST OF TABLES viii 

INTRODUCTION 1 

GEOMAGNETIC FIELD LINES IN SPACE 4 

PARTICLE MOTIONS AND THE UPPER ATMOSPHERE 9 

RADIATION DAMAGE TO SPACECRAFT 14 

RADIATION HAZARDS TO HUMANS 18 

THE FUTURE OF THE SOUTH ATLANTIC ANOMALY 19 

REFERENCES 27 

APPENDIX 28 


The Geomagnetic Field... 


v 




FIGURES 


Figure 1. Contours of total field intensity, in nT, at the Earth’s surface, using IGRF 95 2 

Figure 2. (Upper) Total Field (nT) due to centered, inclined dipole. (Lower) Total Field (nT) due 
to offset, inclined dipole 3 

Figure 3. Contours of total field intensity, in nT, at various altitudes above the Earth's surface, 
using IGRF 95. Upper left 200 km altitude, upper right 500 km, lower left 1000 km, lower right 
5000 km 3 

Figure 4. Field strength as function of altitude above four anomalous areas of the Earth’s surface ... 4 

Figure 5. Lines of force from opposite sides of the Earth. Lines originate from 20 S, 40 S, and 
60 S, and from 40 W (near the SAA) and 140 E (near the Siberian and Antarctic highs) 5 

Figure 6. Lines of force from 70 S, for various longitudes 6 

Figure 7. Lines of force originating from 10 S, 30 S, and 50 S, and from each 30 degrees of longi- 
tude, as projected on the Earth’s surface 7 

Figure 8. Lines of force crossed by the shuttle on an orbit crossing the SAA. On the left is an orbit 
with 28.5 degrees inclination. On the right is an orbit with 57m degrees inclination 8 

Figure 9. (Left) Field Strength vs latitude for starting latitudes 20 S, 30 S, and 40 S, all at 60 W 
longitude. (Right) Field strength vs latitude for starting longitudes 0, 30 W, and 60 W 9 

Figure 10. Heights of lines of force starting at 20 S, 30 S, and 40 S, all along longitude 60 W. The 
gyroradius for electrons (e), and protons (p), are given in meters for 1,10, and 20 Mev particles 
at the equator with pitch angles of 90 degrees 10 

Figure 11. Height of field lines vs latitude for lines starting along longitude 53 W. Numbers on lines 
are mirror points for particles with equatorial pitch angles indicated. Particles are scattered or ab- 
sorbed in the atmosphere below about 100 km altitude 1 1 

Figure 12. Normalized particle flux as function of pitch angle and anisotropy index n, at equator.. 12 

Figure 13. Number density profile with height for atomic oxygen (O), molecular nitrogen (N,), 
and molecular oxygen (CL) 1 3 

Figure 14. Location of Sudden Event Upsets (SEU’s) for TOPEX for the years 1992-1998 with 
contours of total field intensity at 1000 km altitude for 1995 15 

Figure 15. Distribution of Ap-index values for times of TOPEX SEU’s 16 

vii 


The Geomagnetic Field... 



Figure 16. Distribution of dose rates as functions of latitude and of longitude for Skylab 
(Dec. 1973-Jan. 1974) and for Mir (Mar. 1995). Mir’s results were adjusted to Skylab’s to 
compensate for altitude difference (after Badhwar (1997)) 19 

Figure 17. Secular variation of total field (in nT/yr) from IGRF model 20 

Figure 18. Secular variation of total field as observed at longer operating geomagnetic 
observatories around the South Atlantic. The location of the observatories is shown in the 
upper left, with secular change from figure 17 superimposed. A least-squares fit of the 
observed data to a straight line is shown with the slope of that line indicated 22 

Figure 19. Total field contour map (in nT) for the years 2025 (upper), 2050 (center), and 2100 
(lower) extrapolated from IGRF 95 24 

Figure 20. Total field contour map (in nT) for the year 2235, extrapolated from IGRF 95 25 

Figure 21. Location of geomagnetic equator (where dip is zero) for years 1995, 2100, and 2235, 
all based on IGRF 95 26 

TABLES 

Table 1. Some LEO spacecraft 17 

Table 2. Secular variation of the total field at observatories around the South Atlantic 21 


viii 


The Geomagnetic Field... 


INTRODUCTION 


It has been known for more than 40 years that the anomalously weak geomagnetic field in the South 
Atlantic Ocean, near the coast of Brazil, is important to the radiation environment near Earth. Radiation 
belt particles, trapped by lines of force of the geomagnetic field, approach closer to the Earth’s surface 
here than elsewhere. This is an area of intense radiation in near-space and is a subject of extensive re- 
search by space scientists. This radiation is usually described on a heuristic basis from measurements 
made in space and “models” have been developed on this basis. There are also attempts to understand the 
radiation in terms of the laws of physics, i.e., particle motion in the geomagnetic field, and how to effect 
shielding from this radiation. Space mission and instrument designers today rely on both these approaches. 

There is an increasing use of near-Earth orbits for spacecraft to study the Earth, for space telescopes, for 
the manned space station, and for military and commercial applications. All these spacecraft and humans 
in space suffer very significant damaging effects of radiation. Since most of this radiation is related to the 
geomagnetic field, this report will take a brief but closer look at the correspondence between the geomag- 
netic field and radiation and how this situation may become worse in the future. 

Many studies of radiation in near space have treated the geomagnetic field as if it was due to a simple 
dipole. With the current field model we examine details of the field at any latitude, longitude, and altitude, 
out to about 4 Earth radii (25,484 km) where solar time-varying effects become important and the mag- 
netospheric boundary is felt. Many features of the actual geomagnetic field in space have not been illus- 
trated and this report attempts to show that. 

The IGRF spherical harmonic coefficients for 1995 (Barton, C.E. et al., 1996) are used to represent the 
actual geomagnetic field, as of the beginning of the year 1995. This version of the IGRF, as of this writ- 
ing, is the most recently approved reference field and has been extensively studied in many publications. 
This IGRF model also carries the time derivative of the spherical harmonic coefficients so one may calcu- 
late the time variations of the field (although one is cautioned that the accuracy will degenerate if the 
calculation is extended too far from 1995). The three-dimensional complexity of the geomagnetic field is 
such that one should have a computer to display any of the many aspects of the field geometry, or show 
how the degree and order of the IGRF terms influence the field configuration 1 . Here we used the IGRF to 
degree (n) and order (m) 10, although fewer terms were occasionally used to see how significantly they 


'A desktop computer was used to calculate the geomagnetic field components from the International Geomagnetic Reference 
Field (IGRF) coefficients and to produce simple graphics. Computer codes for calculations and graphics were written in 
TrueBasic. A subroutine to calculate the field from the spherical harmonic coefficients follows that of Quinn (Quinn, et al., 
1995) and sources quoted there, converted from the original Fortran to TrueBasic language (see appendix). These programs 
were run on a Power Macintosh 9600 desktop computer. 

There are several programs available to calculate the field components from the IGRF coefficients, for example see http:// 
www.ngdc.noaa.gov/IAGA/wg8/wg8.html. 

To make global maps of the geomagnetic field the total field was calculated for each 5 degrees of latitude and 5 degrees of 
longitude and compared to published tables of these values, where available. The results matched to within a few nT. This 
matrix of values was entered in a graphics program (Spyglass) and the results contoured to produce maps of the total field of 
the Earth. 
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changed the results and to see how the field appears with only the tilted and displaced dipole. 

It is the total field magnitude, rather than field components, that are sufficient for purposes of this report; 
however, the directions of the field are used in the ray tracing of the lines of force. 

Four principal anomalous regions of the field at the Earth’s surface are the localized highs in northern 
Siberia, in northern Canada and near Antarctica south of Australia, and the large, elongated low near the 
coast of Brazil. This later feature is the South Atlantic Anomaly (SAA) 2 , with a minimum of about 23,000 
nT located at 26 S, 53 W, about 700 km inland from the coast of Brazil. It is frequently said that the SAA 
is due to the displacement of the geomagnetic dipole axis from the center of the Earth, in the direction of 
the Pacific Ocean and the tilt of the dipole with respect to the rotational axis. 

This point is made clear in figure 2. The upper part of this figure shows the field that would result from an 
axial, inclined dipole. That is the field that is produced using only the first three terms of the IGRF 1995, 
all of which have n = 1 . In equatorial latitudes the contours follow the magnetic equator. 

The lower part of figure 2 shows the field that would be produced by an offset, inclined dipole. This field 
is due to only the first 8 terms of the IGRF 1995, where terms have n = 1 and 2. The displacement of the 
dipole is 526.89 km in the direction of the northwest Pacific (21.47 N, 144.77 E). The inclined, offset 
dipole creates an Antarctic high south of Australia and an SAA off Brazil but it does not create the highs 
in Northern Canada and Siberia. The SAA in this figure has many of the characteristics of the SAA of the 
complete field (figure 1). However, the actual SAA is stretched more towards the southern tip of Africa, 
showing that terms with n greater than 2 play a significant role shaping that Anomaly. The shape of the 
field contours in figure 2 (Lower) greatly resemble the actual field at an altitude of 5000 km (figure 3, 
lower right), showing that terms with n greater than 2 have little effect at this altitude. 


2 The SAA was also called the Brazilian Anomaly (Dessler and Karplus, I960) or Capetown Anomaly (Heirtzler and Hirshman, 
1960: Cain, etal., 1968). 
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Figure 2. (Upper) Total Field (nT) due to centered, 
inclined dipole. (Lower) Total Field (nT) due to offset, 
inclined dipole. 
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Figure 3. Contours of total field intensity, in nT, at various altitudes above the Earth’s 
surface, using IGRF 95. Upper left 200 km altitude, upper right 500 km, lower left 1000 
km, lower right 5000 km. 
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Maps of the total field at elevations of 200, 500, 1000, and 5000 km are shown in figure 3. Many low- 
Earth orbit (LEO) spacecraft fly at altitudes between 500 and 1000 km and the field they experience is 
shown here. The shape of the contours of the SAA are rather different than at the Earth’s surface. The 
weakest part of the field is still off the coast of Brazil but there is a tendency for the low to be stretched in 
the direction of Africa. From 200 to 500 km the strength of the SAA is reduced by about 4000 nT, and 
from 500 to 1000 km the SAA is reduced by another 4000 nT. At 5000 km the overall field strength is 
much reduced over its value at lower elevations, rarely reaching 10,000 nT, and the difference between 
highs and lows is less pronounced. 



( 1ER) (2ER) (3ER) (4ER) 

Altitude above Earth's Surface (km) 


Figure 4. Field strength as function of 
altitude above four anomalous areas of 
the Earth’s surface. 


The variation of the field with altitude over the four anomalous regions is shown in figure 4. The SAA low 
covers a larger portion of the Earth’s surface than the highs and figure 4 shows that the variation with 
height for the SAA is much less than for the other three areas. For all the regions the field strength at 5000 
km is about one-third what it is on the surface. At 1 0,000 km (somewhat less than 2 Earth radii) the field 
strength is about one-tenth the surface value. Above about 15,000 km (not shown) the field strength is 
much reduced and the four anomalous regions (including the SAA) have lost their distinctive appearances. 
Geostationary orbits (36,000 km) are beyond the influence of the SAA, although there are other radiation 
hazards, due to other effects, at this altitude. 

GEOMAGNETIC FIELD LINES IN SPACE 

Usually early drawings of the Earth’s magnetic field lines were more cartoon than fact and did not reflect 
their true configuration. However, there were attempts to describe the geomagnetic field in a fashion 
useful for studying the radiation environment (for example Mcllwain, 1961; Mlodnosky and Helliwell, 
1962). With the availability of inexpensive but fast desktop computers it is easy to get an accurate presen- 
tation of the field using all the spherical harmonic components. We show some examples here. 
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The straightforward procedure for drawing a line of force (ray tracing) is to start with a point at a given 
latitude, longitude and altitude. For this point the IGRF is used to find the direction of the field (declina- 
tion and inclination). The next point will be a specified distance in the direction indicated by the field. The 
accuracy of the position of the line of force will depend upon this step size. The smaller the step size the 
more accurate the line but the longer it will take the program to run. Using different step sizes we found 
that a 1 10-km step produced sufficient accuracy, although this was slightly degraded along arcs at dis- 
tances of about 4 Earth radii. 

For this second point the direction of the field was again calculated, and position of the third position 
determined, etc., with the process stopping when the next altitude would be less than zero. This provides a 
series of points defining a line of force with the coordinates and field components specified for each point. 
Usually we took lines originating in the Southern Hemisphere (because the SAA is in the Southern Hemi- 
sphere) with a loop stepping the origin to different starting latitudes and/or longitudes as each line is 
finished. In special cases, for example, where we are plotting lines crossing the shuttle’s path (figure 7), it 
was necessary to start the line on the shuttle’s path and proceed in both south and north directions from there. 

To illustrate the asymmetry of the lines of force we show (figure 5) lines of force originating along 40 W 
near the longitude of the SAA and also along 140 E, which is separated by 180 degrees and is near the 
longitude of the Antarctic and Siberian highs. Here the lines originate at 20, 40, and 60 S. Lines which 
originate at 60 S reach nearly 4 ER near the equator along 140 E but only about 1 ER on the opposite side 
of the Earth. Notice that the line which leaves 60 S terminates near 80 N along 140 E, but only at about 45 
N along 40 E. The tilt of the dipole axis is evident. When the Earth rotates during the course of the day, 
this asymmetric field creates a changing, magnetically complex region in near space. 


Lines along 40 W I Lines along 140 E 



Distance East-West (km) 

Figure 5. Lines of force from opposite sides of the Earth. Lines 
originate from 20 S, 40 S, and 60 S, and from 40 W (near the 
SAA) and 140 E (near the Siberian and Antarctic highs). 
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Figure 6 is another display of field lines, originating at various longitudes but always from 70 S. This 
displays the range of altitudes reached by lines originating from different longitudes. For this latitude (70 
S) lines starting from 80 W are closest to the surface of the Earth and the line staring from 120 E are the 
most distant. It is clear that any charged particle within 1 or 2 ER would experience different forces 
depending upon which side of the Earth it approached. 
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Another way of getting a broad view of field lines is to show their projection on the surface of the Earth 
(figure 7). Here lines originate at 10, 30, and 50 S and at various longitudes. The projection of the line 
which originates near the SAA (30 S, 30 W) appears to cross the line originating at 50 S, 30W. Other lines 
appear to cross as well. 



Figure 7. Lines of force originating from 10 S, 30 S, and 50 S, and 
from each 30 degrees of longitude, as projected on the Earth’s surface. 
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Actually they are at different altitudes and do not physically cross in space. If one looks at the lines of 
force crossing the shuttle’s path one gets another perspective (figure 8). Since the shuttle is at constant 
altitude (480 km here) it will cross more polar or more equatorial lines as it moves from high to low 
latitudes. These are shown for two common orbital inclinations: 28.5 and 57 degrees. In both cases an 
orbit was chosen which goes through the SAA. On the 28.5 degree inclination orbit none of the lines of 
force reach high (auroral) latitudes. With the 57 degree orbit the shuttle is frequently near the ends of a 
line of force. In those higher latitudes the lines of force are nearly vertical and do not extend much past the 
flight path. At the equator the shuttle is near the top (maximum altitude) of the line of force and the lines 
of force are more horizontal. 
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Figure 9 is yet another illustration of the lines near the SAA. It shows the field strength along three lines 
of force starting at latitudes 20 S, 30 S, and 40 S, all along 60 W and along three lines originating at 0, 30 
W, and 60 W and all at 30 S. The further south the origin, the weaker the field near the equator because of 
its higher altitude, but the field is stronger at the ends of the lines where the lines are at higher latitudes. 
The field strength at the equator of the line starting on the Greenwich meridian is less than the others 
because that line rises to higher altitudes along the equator. One notices that all lines have the same 
strength at about 28 S, very near the latitude of the SAA. Of course, this does not mean that all lines have 
the same altitude at the crossover point shown. 




Figure 9. (Left) Field strength vs latitude for starting latitudes 20 S, 30 S, and 40 S, all at 60 W 
longitude. (Right) Field strength vs latitude for starting longitudes 0, 30 W, and 60 W. 


PARTICLE MOTIONS AND THE UPPER ATMOSPHERE 


For human safety and for engineering purposes it is important to know the radiation environment expected 
in near-space missions. Ad hoc “models” have been developed from observational data. There are two 
models. One model describes the electron flux and another the proton flux. NASA Trapped Radiation 
Models AE-8 and AP-8, for electron- and proton-produced radiation respectively, are available from the 
National Space Science Data Center (http://nssdc.gsfc.nasa.gov) and described by Fung, 1996. Countries 
other than the U.S., notably Britain and Russia, have other models. The models give the particle flux as a 
function of particle type, particle energy, position, and time. 

No extensive discussion of charged particles in space will be given here, only enough of the basic prin- 
ciples to understand the effect of the geomagnetic field on the radiation. A more complete and convenient 
description can be found on the homepage of the (Belgium) Space Environment Information System 
( http://www.spenvis.oma.be/spenvis/) 
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The motion of a charged particle is characterized by the particle’s energy, mass, charge and pitch angle at 
the equator (angle made with the line of force). Their motion in the radiation belts is depicted as being of 
three basic types, all occurring simultaneously and all controlled by the geomagnetic field (for general 
discussions see Walt, 1994, or Roederer, 1970). The first type of motion is a circular motion with a 
gyroradius about the field line (period of the order of milliseconds). This radius depends upon the strength 
of the field where the particle is located and its local pitch angle. Stronger fields cause the radius to be 
less; the radius is greater for greater pitch angles and greater energies. Figure 10 shows the radius for 
selected points, selected field lines, and selected energies for electrons and protons. It is seen that these 
radii are of the order of meters for electrons and kilometer for protons. The radius will change as the 
particle moves along a field line into different field strengths. 



Figure 10. Heights of lines of force starting at 20 S, 30 S, and 
40 S, all along longitude 60 W. The gyroradius for electrons (e), 
and protons (p), are given in meters for 1,10, and 20 Mev 
particles at the equator with pitch angles of 90 degrees. 


The second motion is the bounce back and forth along a field line from one hemisphere to the other, 
reversing direction at a mirror point. This period is on the order of seconds. This motion is an adiabatic 
invariant, requiring that the magnetic field strength divided by the square of the sine of the pitch angle is a 
constant anywhere along the field line, regardless of the charge, mass, or energy of the particle. Since the 
pitch angle at the mirror point (O m ) is 90 degrees we can find the value of the field strength at the mirror 
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point (B m ) by the equation 


.2 .2 

Bg /sin (Og) — B m /sin (Oj^) = B^ 


( 1 ) 


where B e is the field strength at the equator and d> e is the pitch angle at the equator and O m is the pitch 
angle at the mirror point (90 degrees). 

Since we know the angles of the field along the field line we can find the place where the dip angle 
changes sign, i.e., the equator, and the value of B e there. Equation 1 shows that there is a different value 
of B m for each equatorial pitch angle. Knowing B m we can find the altitude of the mirror point. Figure 1 1 
shows the heights of lines of force vs latitude for south points of origin along longitude 53 W. Only the 
lower 1000 km of these lines are shown and mirror altitudes are shown along each line for various equato- 
rial pitch angles. Although the longitude of these lines goes through the SAA, this figure does not give any 
indication why charged particles are trapped at the SAA. The longer (higher altitude) field lines have 
lower pitch angles near the Earth’s surface. 


Mirror points along 53 W longitude 



Figure 11 . Height of field lines vs latitude for lines starting 
along longitude 53 W. Numbers on lines are mirror points for 
particles with equatorial pitch angles indicated. Particles are 
scattered or absorbed in the atmosphere below about 100 km 
altitude. 
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It should be kept in mind that the particle flux density is not the same for all angles at the equator. Here 
one relies on the cumulative observations that have been made on many spacecraft. The distribution of 
particle flux at the equator is expressed as function of pitch (O) angle by the expression 

J(O) = J sin (<J>) (2) 


Where J is a constant (value of J(O) when 0=90) and n is the anisotropy index and has been given by 
Fung, 1996. Figure 12 is a plot of sin 11 (O), or the ratio of J(d>)/J, vs O for values of n from 1 to 8, indicat- 
ing that, even for relatively small values of n, most of the flux will be at very high pitch angles. 



Fung shows the electron and proton particle density of a function of n, but the number of particles with 
lower values of n is relatively small. Then, referring to figure 1 1 we see that most particles have a mirror 
location in the upper part of their field lines and much fewer come to lower elevations where they might 
interact with the atmosphere. 

The third type of motion is the drift of particles around the Earth from east to west or west to east, usually 
taking about an hour. As particles gyrate about a field line they reach neighboring field lines to the east or 
west. There may be a stronger or weaker field on the neighboring line, shortening or lengthening the 
gyroradius. This will cause the particle to drift to the east (electrons) or west (protons). This effect may be 
different at different altitudes, latitudes, and longitudes because the field strength varies with these vari- 
ables. This motion of charged particles is called the “ring current.” 
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Scattering by the atmosphere removes particles from the lower end of the radiation belts. This is illustrated 
by the reduction in astronauts’ dose rate during sunspot maximum, compared to sunspot minimum 
(Badhwar, et al., 1997). At sunspot maximum the atmosphere is believed to expand due to heating, remov- 
ing more particles from the lower end of the radiation belts. Although there are more solar flares — produc- 
ing more particles in the radiation belts — during solar max, this effect is of shorter duration and is counter 
to the effect of removing radiation-belt particles by the expanded atmosphere. 

Models of the upper atmosphere are available on-line from the National Space Science Data Center 
(NSSDC) at URL ( http://nssdc.gsfc.nasa.gov/space/model/atmos/msice.html ). This source provides data 
for the MSIS-E-90 (Mass-Spectrometer-Incoherent-Scatter) model, which in turn is based largely on the 
COSPAR International Reference Atmosphere (CIRA) 1986 model. This model is an improvement over 
the frequently quoted U.S. Standard Atmosphere (1976). Figure 13 shows, as an example, a height-density 
profile of the atmosphere for 06 hrs UT, 1 January 1995, in the vicinity of the SAA. 



On a worldwide basis the atmospheric density is highest in October and April and lower in July and 
January. It also varies with latitude, longitude, time of day, geomagnetic activity and, as just mentioned, 
sunspot cycle. NSSDC atmospheric models, plotted for different locations and times, show that most of 
the density change occurs in the lower atmosphere and that figure 13 is a typical example. 

The number density is dominated by molecular nitrogen below 80 to 100 km. Above this altitude atomic 
oxygen dominates the particle number density. However, the total particle density (O plus NA drops off 
drastically by about 7 orders of magnitude. The sensible top of the atmosphere is usually taken to be 100 km. 
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Interactions with the atmosphere are not necessarily responsible for the increase in radiation over the 
SAA. However, the atmosphere is responsible for removing particles from the radiation belts. 

RADIATION DAMAGE TO SPACECRAFT 

There are many types of hazards to spacecraft. Aside from the hazards caused by space debris and the 
problems associated with the vacuum of space, there are effects caused by the interaction of spacecraft 
with its plasma environment. These include spacecraft wake effects, spacecraft contamination effects, 
spacecraft power system effects, and spacecraft electromagnetic effects. There are major effects due to the 
natural, undisturbed radiation in space and these we explore here. 

All spacecraft suffer radiation damage. Designing electronics and other equipment to minimize the effect 
of this damage is a major engineering concern, affecting both the choice and redundance of electronic 
components and circuits, as well as software design. Especially over the SAA, there is significant radia- 
tion damage to electronic, optical, and computer systems. This radiation is of various types, depending 
upon the type of particle and its energy, the amount of radiation shielding, and the sensitivity of the mate- 
rials receiving this radiation. One common type of radiation damage is spacecraft charging whereby 
charge is built up with time until arcing occurs. Another and very common type of damage is a “Single 
Event Upset” or SEU, whereby a single and critical digital component changes state. Clever design can 
frequently cause such problems to be reduced, although the spacecraft may be inoperative for a period of 
time. Some spacecraft turn off electrical circuits or optical systems when passing over the SAA. 

Most ground controllers keep records of the type of damage, the time and where such damage occurred. 
An archive of such records is available from the National Geophysical Data Center in Boulder, URL 
address ( http://spider.ngdc. noaa.gov:800/production/htmI/GQES/anom5jt.txt ) or by ftp :// 
ftp.ngdc.noaa.gov/lSTP/ANOMALIES/. 

An example of the SEU’s for one spacecraft is given in figure 14. This is for TOPEX, at an altitude of 
1340 km, for the period 1992-1998, with total field strength contours shown for approximately the same 
altitude. The points do not cluster evenly within a contour but are more scattered in a northeast-southwest 
direction. (Characteristically, SEU’s for various spacecraft are not rigidly confined by any total field 
contour line.) 
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There were 282 SEU’s over this 6-year period. This rate of a few per week to a few per month seems 
typical. It is seen that most of the SEU’s are located in and about the SAA, although there is no clear 
contour that defines where they will occur. An examination was made of South America geomagnetic 
observatory records for times when there was an SEU nearby. These records indicated no unusual geo- 
magnetic activity. Another investigation looked at the 3-hour Ap activity index when there was an SEU. 
Figure 15 shows the distribution of Ap when there was an SEU. There are relatively few high Ap and the 
distribution of Ap is not greatly different from the distribution of all Ap during the 1992-1998 period. 
Although it is known that there is a burst of radiation in space when there are solar flares it is impossible 
to predict exactly when an SEU will occur. Table 1 is a list of some spacecraft which have been, are, or 
will be, in Low-Earth Orbit and suffer radiation damage. A notable case is the Hubble Space Telescope, 
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whose instruments are turned off over the SAA to prevent damage to its optical fibers, electronic compo- 
nents, and solar arrays. Other spacecraft have had problems with star cameras and other instruments. The 
shuttle, Mir, and the International Space Station Alpha are manned space missions which have operated or 
will operate for extended periods of time in this environment. Spacecraft in higher orbits also suffer 
radiation damage. For example, spacecraft in geosynchronous orbits, although at an altitude of 35,790 km 
(about 6 Earth radii) and in an equatorial orbit also suffer radiation damage. These spacecraft rotate with 
the line of force on which they are located and these lines terminate at very high latitudes on the Earth’s 
surface. 



Ap value 

Figure 15. Distribution of Ap-index values for 
times of TOPEX SEU’s. 
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Table 1 . Some LEO Spacecraft. 


Name 

Country 

EXISTING 


ASTRID-2 

SWEDEN 

COBE 

US 

ERS-1 

ESA 

EUVE 

US 

GRO 

US 

HST 

US 

LANDSAT-5 

US 

LANDSAT-7 

US 

MIR 

RUSSIA 

NOAA-12 

US 

OERSTED 

DEN. 

ROSAT 

GER. 

SAMPEX 

US 

SEAWTFS 

US 

SHUTTLE 

US 

SPOT 

FRANCE 

SUNSAT 

S. Afr. 

TOPEX 

US/FR 

TRIMM | 

US 

UARS 

US 

PLANNED 


CHAMP 

GER. 

DMSP-B5 

US 

FEDSAT 

AUST. 

GRADSAT 

DEN. 

NPOESS 

US 

SAC-C 

ARGEN. 

SAC1-1 

BRAZIL 

VCL 

US 

GLAS 

US 


Launch 

Inch 

Mean 

Date 

(deg.) 

Alt. (km) 

12/98 

83 

1000 

11/89 

99 

882 

7/91 

99 

784 

7/92 

28 

506 

4/91 

28 

434 

4/90 

28 

595 

4/84 

98 

703 

4/99 

98 

700 

2/86 

52 

392 

5/90 

99 

815 

2/99 

97 

715 

6/90 

53 

544 

7/92 

82 

593 

9/97 

98 

705 

many 

28.5, 57 

300 

several 

99 

825 

2/99 

96 

715 

6/92 

66 

1338 

11/97 

35 

350 

8/91 

57 

578 

/OO 

85 

375 

8/99 

polar 

840 

/01 

70 

750 

/ 02 

polar 

600 

/08 

polar 

833 

12/ 99 

98 

702 

7/99 

99 

750 

2/00 

65 

400 

/01 

94 

600 
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RADIATION HAZARDS TO HUMANS 


Radiation in the radiation belts offer several constraints and hazards to humans in low-Earth orbit. The 
ionizing radiation found in LEO is thought to come, largely, from two sources: (1) trapped radiation which 
originates on the Sun and with the proton component being the most dangerous, and (2) galactic cosmic 
rays, which are extremely energetic and which occur at unpredictable times. Shielding offers little protec- 
tion from radiation belt particles from the Sun and virtually no protection from galactic rays (see Wilson, 
et al., 1997). Outside the Earth’s magnetic field (outside the magnetosphere) there are other very signifi- 
cant radiation hazards but, since those are not controlled by the Earth’s magnetic field, they will not be 
discussed here. The environment of low-Earth orbit is where humans will be spending much time in the 
coming years and decades. Within NASA the Spacecraft Radiation Health Program at Johnson Space 
Center t http://sn-io.jsc.nasa.gov/sreg/techmemo.htm') leads in the research and development activities that 
address the health effects of space radiation exposure to astronauts. Astronauts carry personal radiation 
dosage badges and radiation monitors are located throughout the shuttle. There is a space radiation analy- 
sis group (SRAG) which predicts, in real time, the exposure that might occur during a solar flare and 
terminates extravehicular activity if that is necessary. 

It is impossible to predict how much radiation an astronaut will receive, although one can make an esti- 
mate of the probability of receiving a given amount. However, there is a chance that these limits will be 
exceeded in space. NASA regulations for limits of exposure, measured in REM’s, call for 15 REM for 30 
day exposure, 50 REM for annual exposure, 1 50 to 400 REM for exposure for men (more REM with 
increasing age), 100 to 300 for exposure for women (more REM with increasing age). 

The radiation environment experienced by astronauts on Skylab and cosmonauts on Mir, as a function of 
their geographic location, has been studied by Badhwar (1997). Figure 16 shows figures from that paper, 
clearly illustrating that there is a major peak in radiation at the SAA, although there are minor peaks in 
other places. The 300 nGy/min peak over the SAA corresponds to about 30 mRad/min (30 mREM/min). 
Figure 8 shows the line of force cut by the shuttle, which is at 480 km altitude, orbital inclination 28.5 
degrees. Although this line is steeply inclined to the vertical, the shuttle is well beyond the sensible atmo- 
sphere. The average value of radiation in the SAA is less than this peak value and is about half of this or 
15 mREM/min. An astronaut would have to spend 1000 minutes (17 hrs) in the SAA during 30 days to 
reach his/her exposure limit. For a U.S. astronaut or Russian cosmonaut who stays in space for a year, 
perhaps on the space station, there is the 50 REM annual limit (note that is much less than 12 times the 
monthly limit). The radiation dosage of astronauts on the space station will be continuously monitored 
from Earth. During their flights to and from the moon the Apollo astronauts were lucky because they 
would have received a lethal dose of radiation if a strong solar event had occurred while they were in 
space — and such disturbances occur several times per year. A recent National Academy of Science study 
(National Academy of Science, 1996) showed that manned flights to Mars, because they are of long 
duration and outside the shielding of the magnetosphere, will not be possible, regardless of wishful think- 
ing, until the shielding problem can be solved. 

It is interesting to note that TOPEX, at an altitude of 1334 km, inclination 63 degrees, spent 42 minutes in 
the SAA over a 3-year period ( 14 min/yr). It received 2.3 rad per orbit or an accumulated dose of 12,000 
rad/yr behind 1 gm/cm 2 Al shielding. This is far in excess of the lifetime limits for humans. 
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Longitude (degrees) L*m«fc (degtee.) 

Figure 16. Distribution of dose rates as functions of latitude and of longitude for Skylab 
(Dec. 1973-Jan. 1974) and for Mir (Mar. 1995). Mir’s results were adjusted to Skylab's to 
compensate for altitude difference (after Badhwar (1997)). 


Shuttle flights with a 28.5 degree inclination do not pass through the entire SAA. Shuttle flights with 57 degree 
inclination do pass through the entire SAA but receive less radiation because they spend less time in the SAA. 

THE FUTURE OF THE SOUTH ATLANTIC ANOMALY 

Since the geomagnetic field is not static it is worthwhile investigating what the SAA and associated 
radiation hazard might be in the future. Figure 17 shows the secular change in total intensity of the Field 
for the year 1995, as given by the IGRF. It shows that the greatest rate of decrease of the field is in the 
South Atlantic southwest of Capetown where the secular variation is more than 100 nT/yr. There is a 
secondary area of rapidly decreasing field in the western North Atlantic where it exceeds 80 nT/yr. Mea- 
surements over the past century, for much of the world, show essentially the same rate of change as for 
1995 (Sabaka, et al., 1997). It is understood that data taken at some specific locations will differ from 
those shown by the spherical harmonic contours since the process of generating the spherical harmonic 
coefficients will possibly smooth over all values to fit a spherical harmonic expression. 
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Precisely extrapolating the rate of change of the geomagnetic field into the future is impossible because 
the reason for change is unknown. Here we can only draw upon what has happened in the past and assume 
that will continue for the next few decades. For some observatories, especially in the Northern Hemi- 
sphere, and in India, the slope of field strength vs time has shown a complete change of slope in the past. 
Many investigators have not tried to predict the field in the face of this erratic behavior. For the SAA 
region we can look at the records from 1 1 observatories that exist around the South Atlantic region. Those 
with records extending for more than 25 years are listed in table 2, and their locations shown in figure 18. 
Table 2 also gives the period of time for which there are yearly mean records and the rate of change of 
total field that those records show, calculated from a least-squares fit to a straight line. These data were 
obtained from NGDC. Figure 18 also shows the data for five observatories with the longer records and 
illustrates that the field is changing at a constant rate. Comparing these rates of change to figure 17 it is 
seen that the change at Hermanus (HER) and Syowa Base (SYO) are greater than figure 17 shows, but 
other observatories are about the same as the global map of figure 17 indicates. 
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Table 2. Secular variation of the total field at observatories around the South Atlantic. 


Annual Means From NGDC Annual Means Web Database: 

Code 

Sec, Var. 
(nT/yr) 

Years 

Name 

HER 

-110.4 

1940-1999 

Hermanus, S. Africa 

KOU 

-105.4 

1958-1975 

Surinam 

LAS 

-46.7 

1965-1995 

Las Acacias, Argentina 

LQA 

-47.4 

1928-1987 

La Quiaica, Argentina 

MBO 

-14.5 

1952-1988 

M’Bour, Senegal 

PIL 

- 57.6 

1905-1997 

Pilar, Argentina 

SYO 

-101.3 

1955-1997 

Syowa Base, Antarctica 

TRW 

- 77.9 

1958-1995 

Trelew, Argentina 

TSU 

-57.7 

1952-1997 

Tsumba, Namibia 

TTB 

-69.2 

1960-1992 

Tatuoca, Brazil 

VSS 

-25.2 

1920-1997 

Vassouras, Brazil 

Annual Means From IGPP Monthly Means Web Database: 

HER 

-87.4 

1992-1994.8 

Hermanus, S. Africa 

HER 

-81.4 

1995-1998.5 

Hermanus, S. Africa 

LAS 

-65.1 

1992-1994.6 

Las Quiaica, Argentina 

TSU 

-59.6 

1992.2-1994.9 

Tsumba, Namibia 

TSU 

-65.6 

1995-1998.5 

Tsumba, Namibia 
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Figure 18. Secular variation of total field as observed at longer operating geomag- 
netic observatories around the South Atlantic. The location of the observatories is 
shown in the upper left, with secular change from figure 17 superimposed. A least- 
squares fit of the observed data to a straight line is shown with the slope of that 
line indicated. 


Also in table 2 we give rates of change calculated from monthly mean values, for some observatories, with 
the data obtained from IPGP Paris ( http://obsmag.ipgp.jussieu.fr/AM-rqobsmag.html '). The monthly mean 
data only extends for the last 8 years and does not follow so smooth a curve as the yearly means. The 
slope of these curves was also determined by a least-square fit to a straight line and those slopes are given 
in table 2. It is seen that the monthly means for Hermanus give values closer to the values in figure 17, 
than does the yearly mean. This suggests that, in the last few years, the field at Hermanus and Syowa Base 
may have been decreasing faster than at earlier times. 

For periods of time longer than shown in table 2 one might think of the geologic record of the paleomag- 
netic field. Kreer et al. (1983) reported on a paleomagnetic study of Argentine lake sediments from each 
of three adjacent lakes located about 600 km northwest of the Trelew geomagnetic observatory. The 
sediments from one of the lakes extended back 14,000 years, but sediments from the other lakes did not 
extend back so far. The declination and inclination records from the three lakes could be stacked to pro- 
duce a record of the direction of the field for the last 6,000 years. The probable error of these directions is 
high although they show a large excursion to the west and up about 500 years ago. 
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The authors show intensity and susceptibility profiles from eight cores from the three lakes but were 
unable to make a stacked intensity profile because of the high variability from lake to lake. From the 
figures shown the intensity seems to have varied erratically by 10 to 20 percent with numerous short 
excursions of intensity, which can easily be artifacts. The authors point out that these are the first paleo- 
magnetic measurements from South America and much remains to be done before measurements like this 
have high reliability. A search of the literature did not reveal any paleomagnetic measurements from lakes 
from South Africa. 

The Deep Sea Drilling Project and the International Program of Ocean Drilling have drilled a number of 
holes in the South Atlantic Ocean floor. Magnetic measurements are usually made on the deep-sea core 
material but, contrary to cores from the North Atlantic, no comprehensive report of paleomagnetic intensi- 
ties has been made for the South Atlantic. In any case, the sedimentation rate in the oceans is usually 
much slower than in lakes, and time scale so compressed, they are not likely to yield information relevant 
to the last few hundred years. 


It is clear that the spherical harmonic model of the field and its secular change smoothes the field and 
cannot be expected to give exactly correct values at all points on the Earth. For purposes of this report we 
will assume that the rate of change as illustrated in figure 17 is correct enough to be used to extrapolate 
tens to a few hundred years in the future. The calculation is carried out by updating each of the spherical 
harmonic coefficients, according to their first-time derivative, calculating the X, Y, and Z components of 
the field and then calculating B from X, Y, and Z. 

Figure 19 shows the total field for the year 1995 extrapolated to the years 2025, 2050, and 2100. The 
principle feature of the geomagnetic field is the change in the SAA. Each few decades it gets larger and 
the field within it gets weaker. A secondary low shows southwest of Capetown. 

One may define the extent of the SAA by the areal extent of the high dose rates for shuttle and Mir (figure 
16). The high dose rates extend from approximately 15 S to 40 S and from 20 W to 70 W at 480 km 
altitude. This corresponds roughly to the 25,000 nT contour on the Earth’s surface. Using this contour to 
identify the SAA, in the year 2100 the SAA will cover most of South America, the southern part of Africa 
and the South Atlantic Ocean south of 25 S to the Scotia Sea and Antarctica. The size of the SAA will 
have increased by a factor of about 4. This suggests that the radiation hazard to humans in space may be 
increased by a corresponding amount. This will mean that astronauts can spend correspondingly less time 
in space. 

In 1995, the focus of the SAA was clearly over the coast of Brazil. In 2100 there will be two foci: one in 
northern Brazil and a second and deeper one southwest of Capetown in the South Atlantic Ocean. 

Radiation levels with this new SAA will be different from the present for two reasons: (1) because the 
field is less and the mirror points of particles lower with more being lost to the atmosphere (greater 
bounce cone loss), and (2) as mentioned before, the area of this low field will be greater than it is now. 

The first effect would probably reduce the level of radiation, while the second effect would tend to distrib- 
ute the radiation over a larger area with the integrated effect increasing the radiation hazard. It is tacitly 
assumed that at present the radiation belts are in approximate equilibrium with the input of particles from 
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galactic and solar sources balancing the number of particles lost by atmospheric scattering. In the future 
the sources would be the same as today while the losses would be greater. This could seriously deplete the 
particle population of the radiation belts. Able and Thome, 1999, have made progress in quantifying the 
factors which create a radiation hazard over the SAA in 1995. A similar analysis for a possible future field 
is needed with attention given to the balance between the source supply and the atmospheric loss. 



It is interesting to continue to extrapolate the field 240 years in the future, to the year 2235. The focus in 
northern South America intensifies but the field strength approaches zero in the South Atlantic (figure 20). 
Naturally, all components of the field (X, Y, H, and Z) must approach zero if the total field approaches 
zero. This does not constitute a global reversal of the geomagnetic field but is a phenomenon confined to 
the South Atlantic, at that time. 

Using the 1995 IGRF secular change in inclination one may calculate how the geomagnetic equator 
(where the field is horizontal) will change with time. Figure 21 shows that the geomagnetic equator in the 
Atlantic will move north by 10 to 15 degrees in the next 100 years and by 15 to 30 degrees in the next 240 
years. In other parts of the world the change is not great, although it may change long-range high-fre- 
quency radio transmissions. 
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APPENDIX 


This computer code, in True Basic, will plot the field strength B or altitude h vs latitude, with mirror 
points indicated (similar to figure 1 1) or it can be modified to make a number of the other plots shown in 
this report. Note that most of this code, after statement 2400, is the subroutine GEOMAGJ to calculate the 
field components for a given latitude, longitude, and altitude. The file “IGRF95” is called for in the sub- 
routine and a sample listing of values in this file is given following the listing of the main program. Also 
note that a large number of lines is to calculate the location of the mirror points and plot them. This part 
can be skipped by inserting a “Go to” statement and this will materially reduce the number of lines. 


10 [Computer code plot B or h vs lat with mrror pts. 

20 [Written in TrueBasic, Nov 1998, JRH 
30 [stores graphics on clipboard 
40 set zonewidth 6 
50 OPTION NOLET 
60 OPTION ANGLE degrees 

70 DIM h(200),lt(200),ln(200),b(200),dp(200),dc(200) [for pts along line of force 
80 DIM bm(100) [for field values at mirror points 

90 DIM hms(100) / ltms(100),lnms(100) / bms(100) [for variables at mirror pts in south hemisphere 
100 DIM hmn(100),ltmn(100), lnmn(100),bmn(100) [for variables at mirror pts in north hemisphere 
110 !OPEN#4: NAME OUTFILE$,CREATE NEW,ORG TEXT 
120 ! ***** initialize for plotting ***** 

130 LIBRARY "PictLib*" 

140 LIBRARY "MacTools*" 

150 call Copy_clipboard 

160 open#2: screen .05,-4,. 05, .6 !to make "trim" fit 

170 xmin=-55 

180 xmax=55 

190 ymin=0 

200 ymax=10000 

210 set window xmin,xmax,ymin,ymax 
220 ! ***** d raw box ***** 

230 plot lines: xmin,ymin;xmax,ymin 
240 plot lines: xmax,ymin;xmax,ymax 
250 plot lines: xmax,ymax;xmin,ymax 
260 plot lines: xmin,ymax;xmin,ymin 
270 ! ***** plot tics on y axis ***** 

280 for n=ymin to ymax step 1000 
290 plot lines: xmin,n;xmin+l,n 
300 plot lines: xmax-l,n;xmax,n 
310 next n 

320 [ ***** plot tics on x axis ***** 

330 for n=xmin to xmax step 5 
340 plot lines: n,ymin,n,ymin+100 
350 plot lines: n,ymax-100,*n,ymax 
360 next n 

370 [ ***** calculate and plot field line ***** 

380 5 
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390 dlat=l !lat increment along path, deg 

400 for flonk=0 to 0 Istep 30 Ithis steps longitude of origin 

410 for nk=-50 to -30 step 10 !this steps latitude of origin 

420 flatmin=nk {starting lat, deg 

430 flatmax=-nk {final lat, deg 

440 altl=0 {starting alt, km 

450 flatl=flatmin 

460 flonl=flonk 

470 k=l 

480 flat— flatl 

490 alt=altl 

500 flon=flonl 

510 ! 

520 ! ***** stepping along path begins here ***** 

530 ! 

540 CALL GEOMAGJ (alt,flat,flon,t,dip,dec) !get value of t,dip, dec 

550 h(k)=alt 

560 lt(k)=flat 

570 ln(k)=flon 

580 b(k)=t 

590 dp(k)=dip 

600 dc(k)=dec 

610 k=k+l 

620 dalt=-110*dlat*tan(dip) 

630 alt=alt+dalt 
640 alt2=alt 

650 if alt2<0 then go to 1050 ! if alt goes into ground we have gone too far 

660 flat2=flatl+dlat 

670 dlon=110*dlat*tan(dec/60) 

680 flon2=flonl+dlon 
690 go to 1000 

700 if flonl>180 then ! change from flon,flat to x,y 
710 xl=flon 1-360 
720 go to 810 
730 else 

740 if flonl<-180 then 

750 xl=flonl+360 

760 go to 810 

770 else 

780 xl=flonl 

790 end if 

800 end if 

810 yl=flatl 

820 if flon2>180 then 

830 x2=flon2-360 

840 go to 930 

850 else 

860 if flon2<-180 then 
870 x2=flon2+360 
880 go to 930 
890 else 
900 x2=flon2 
910 end if 
920 end if 
930 y2=flat2 
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940 if x2>xmax or x2<xmin then go to 1000 lout of bounds 
950 if y2>ymax or y2<ymin then go to 1000 lout of bounds 
960 base=abs(xl)+abs(x2) 

970 if base>=358 then go to 1000 !to prevent line from jumping from one side to the other 

980 !if alt2<600 then box circle xl-.5,xl+*5,yl-.5,yl+.5 lemphase line is alt is low 

990 plot lines: xl,yl;x2,y2 Iplot line segment 

1000 flatl=flat2 ladvance points 

1010 flonl=flon2 

1020 altl=alt2 

1030 go to 480 

1040 ! ***** b(n) or h(n) vs lt(n) from n=l to kmax ***** 

1050 kmax=k-l 
1060 xxl=-lt(l) 

1070 yyl=h(l) 

1080 for n=2 to kmax 
1090 xx2=-lt(n) 

1100 yy2=h(n) 

1110 plotlines: xxl,yyl;xx2,yy2 
1120 xxl=xx2 
1130 yyl=yy2 
1140 nextn 

1150 ! ***** f ind values for dip equator ***** 

1160 ! 

1170 for k=l to kmax-1 Ithis loop ends on line 1470 
1180 if dp(k)<=0 and dp(k+l)>0 then 
1190 factor=abs(dp(k)/ (dp(k+l)-dp(k))) 

1200 if h(k+l)<=h(k) then 

1210 heq=h(k)-factor*(abs(h(k+l)-h(k))) lalt decreasing or no change 
1220 else 

1230 heq=h(k)+factor*abs((h(k+l)-h(k))) lalt increasing 

1240 end if 

1250 lteq=lt(k)+factor*(lt(k+l)-lt(k)) Hat is always increasing 
1260 if ln(k+l)<=ln(k) then 

1270 lneq=ln(k)-factor sf abs((ln(k+l)-ln(k))) llong is decreasing or no change 
1280 else 

1290 lneq=ln(k)+factor*abs((ln(k+ l)-ln(k))) llong is increasing 

1300 end if 

1310 if b(k+l)<=b(k) then 

1320 beq=b(k)-factor*abs((b(k+l)-b(k))) lb is decreasing or no change 
1330 else 

1340 beq=b(k)+factor*abs((b(k+l)-b(k))) !b is increasing 

1350 end if 

1360 if dp(k+l)<=dp(k) then 

1370 dpeq=dp(k)-factor*abs((dp(k+l)-dp(k))) Idip is decreasing or no change 
1380 else 

1390 dpeq=dp(k)+factor*abs((dp(k+l)-dp(k))) Idip is increasing 

1400 end if 

1410 if dc(k+l)<=dc(k) then 

1420 dceq=dc(k)-factor*abs((dc(k+l)-dc(k))) [declination is decreasing or no change 
1430 else 

1440 dceq=dc(k)+factor*abs((dc(k+l)-dc(k))) I declination is increasing 

1450 end if 

1460 end if 

1470 next k 

1480 ! 
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***** C al C south and north mirror pts for eq pitch angle alpha ***** 


1490 ! 

1500 ! 

1510 for alpha=80 to 30 step -10 [this loop ends at statement 2140 
1520 bm(aipha)=beq/(sin(alpha)*sin(alpha)) lvalue of bm for this alpha 
1530 ! 

1540 ! ***** search for and plot mirror pts, first south of eq, then north ***** 

1550 ! 

1560 for k=l to kmax-1 [start following along line of force, ends at statement 2130 
if dp(k)<=0 then [south of mag equator 

1580 if b(k)>=bm(alpha) and b(k+l)<bm(alpha) then lvalue found, calcuate 
1 590 factor=abs( (bm(alpha)-b(k)) / (b(k+ 1 )-b(k)) ) 

1600 if h(k+l)<=h(k) then 

1610 hms(alpha)=h(k)-factor*(abs(h(k-i-l)-h(k))) !h is decreasing 

1620 else 

1630 hms(alpha)=h(k)+factor*abs((h(k+l)-h(k))) !h is increasing 

1640 end if 

1650 ltms(alpha)=lt(k)+factor*(lt(k-i-l)-lt(k)) 

1660 if ln(k+l)<=ln(k) then 

1670 lnms(alpha)=ln(k)-factor*abs((ln(k+l)-ln(k))) !ln is decreasing 

1680 else 

1690 lnms(alpha)=ln(k)+factor*abs((ln(k+l)-ln(k))) tin is increasing 

1700 end if 

1710 if b(k+l)<=b(k) then 

1720 bms(alpha)=b(k)-factor*abs((b(k+l)-b(k))) !b is decreasing 

1730 else 

1740 bms(alpha)=b(k)+factor*abs((b(k+l)-b(k))) !b is increasing 

1750 end if 

1760 xl=ltms(alpha)-l [plot lat, h box around south mirror pt 

1770 x2=ltms(alpha)+l 

1 780 y 1 =hms(alpha)-50 

1790 y2=hms(alpha)+50 

1800 if hms(alpha)<>0 then box circle -x2,-xl,yl,y2 [draw box 1 

1810 else lvalue not found 

1820 go to 2130 [continue search south of mag equator 

1830 end if lend calc for south mirror pt 
1840 else [north of mag equator 

1850 if b(k)<bm(alpha) and b(k+l)>=bm(alpha) then lvalue found, calculate 

1860 factor=abs((bm(alpha)-b(k))/(b(k+l)-b(k))) 

1870 if h(k+l)<=h(k) then 

1880 hmn(alpha)=h(k)-factor*(abs(h(k+l)“h(k))) !h is decreasing 

1890 else 

1900 hmn(alpha)=h(k)+factor*abs((h(k+l)-h(k))) !h is increasing 

1910 end if 

1920 ltmn(alpha)=lt(k)+factor*(lt(k+l)-lt(k)) 

1930 if ln(k+l)<=ln(k) then 

1940 lnmn(cdpha)=ln(k)-factor*abs((ln(k-i-l)-ln(k))) !ln is decreasing 

1950 else 

1960 lnmn(alpha)=ln(k)+factor*abs((ln(k+ l)-ln(k))) !ln is increasing 

1970 end if 

1980 if b(k+l)<=b(k) then 

1990 bmn(alpha)=b(k)-factor*abs((b(k+l)-b(k))) !b is decreasing 

2000 else 

2010 bmn(alpha)=b(k)+factor*abs((h(k+l)-b(k))) lb is increasing 

2020 end if 

2030 xl=ltmn(alpha)-l [plot lat, h box around north mirror pt 
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2040 x2=ltmn(alpha)+l 

2050 y 1 =hmn(alpha)-50 

2060 y2=hmn(alpha)+50 

2070 if hmn(alpha)<>0 then box circle -x2,-xl,yl,y2 !draw box 2 

2080 else lvalue not found, continue north of mag equator 
2090 go to 2130 

2100 end if lends calc north mirror pt 

2110 go to 2140 ljump out of this loop and finish this search 

2120 end if lends test of dip 

2130 next k lends checking line of force for this value of alpha 

2140 next alpha 

2150 for alpha= 80 to 40 step -10 Iplot marks for mirror pts 
2160 if hms(alpha)<>0 then 

2170 xxl=-ltms(alpha)-.5 Iplot lat, b box around south mirror pt 

2180 xx2=-ltms(alpha)+.5 

2190 yyl=bms(alpha)-300 

2200 yy2=bms(alpha)+300 

2210 Ibox circle xxl,xx2,yyl,yy2 Idraw box 3 

2220 else 

2230 go to 2250 

2240 end if 

2250 if hmn(alpha)<>0 then 

2260 xx3=-ltmn(alpha)-.5 Iplot lat, b box around north mirror pt 

2270 xx4=-ltmn(alpha)+.5 

2280 yy3=bmn(alpha)-300 

2290 yy4=brrtn(alpha)+300 

2300 Ibox circle xx3,xx4,yy3,yy4 Idraw box 4 

2310 else 

2320 go to 2340 

2330 end if 

2340 next alpha 

2350 nextnk 

2360 next flonk 

2380 call Copy_done 

2390 close #2 

2400 END 

2410 ! ***** END OF PROGRAM ***** 

24201 

2430 1 ***** subroutine to calculate field components ***** 

2440 SUB GEOMAG](alt,flat,flon,t,dip,dec) 

2450 ! subroutine GEOMAGJ to calc geomag field components from coefs 

2460 Iread in coefs and sec change, after sub getshc of NGDC 

2470 option nolet 

2480 option angle degrees 

2490 dim gh(224) 

2500 dim dgh(224) 

2510 ddate=2100.0 

2520 infile$="IGRF95" 

2530 factor=ddate-1995.0 

2540 nmax=10 

2550 open #3:name infile$ 

2560 i = 0 

2570 for nn = 1 to nmax 
2580 for mm = 0 to nn 

2590 input#3:m,n,g,h,dg,dh 
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2600 

i = i + 1 

2610 

gh(i) = g 

2620 

dgh(i)= dg 

2630 

if (m <> 0) then 

2640 

i = i + 1 

2650 

gh(i) = h 

2660 

dgh(i)=dh 

2670 

end if 

2680 

next mm 

2690 

next nn 

2700 

close #3 


2710 
2720 
2730 
2740 
2750 
2760 
2770 
2780 
2790 
2800 
2810 
2820 
2830 
2840 
2850 
2860 
2870 
2880 
2890 
2900 
2910 
2920 
2930 
2940 
2950 
2960 
2970 
2980 

2990 dim sl(14), cl(14) ! here set for nmax=14 
3000 dim p(119), q(119) 

3010 dim ext(3) 

3030 !igdgc=2 [geocentric coord., elevkm=dist from earth center 
3040 igdgc=l ! geodetic coord., elevkm=elev above mean sea level 
3050 elevkm-alt 
3060 erad=6371.2 

3070 a2=40680925. 

3080 b2=40408588. 

3090 dtr = .01745329 
3100 r = elevkm 
3110 slat = sin (flat) 

3120 if (90. - flat) < .001 then 

3130 aa = 89.999 

3140 1 300 ft from n. pole 


input: 

igdgc - indicates coordinate system used; set equal to 
1 if geodetic, 2 if geocentric 
flat - north latitude, in degrees 
flon - east longitude, in degrees 

elevkm - elevation above mean sea level (igdgc=l), or 
radial distance from earth's center (igdgc=2) 
erad - value of earth's radius associated with the sh 
coefficients, in same units as elevkm 
a2,b2 - squares of semi-major and semi-minor axes of 
the reference spheroid used for transforming 
between geodetic and geocentric coordinates or 
components 

nmax - maximum degree and order of coefficients 
gh - schmidt quasi-normal internal spherical 
harmonic coefficients 

output: 

x - northward component 

y - eastward component 

z - vertically-downward component 

based on subroutine 'igrf' by d. r. barraclough and 
s. r. c. malin, report no. 71/1, institute of geological 
sciences, u.k. 
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3150 else if (90. + flat) < .001 then 

3160 aa = -89.999 

3170 ! 300 ft from s. pole 

3180 else 

3190 aa = flat 

3200 end if 

3210 clat = cos (aa) 

3220 sl(l) = sin (flon) 

3230 cl(l) = cos (flon) 

3240 x = 0. 

3250 dx=0. 

3260 y = 0. 

3270 dy=0. 

3280 z — 0. 

3290 dz=0. 

3300 sd = 0. 

3310 cd = 1. 

3320 n = 0 
3330 1 = 1 

3340 m = 1 

3350 npq - (nmax * (nmax + 3)) / 2 

3360 if igdgc = 1 then 

3370 aa = a2 * clat * clat 

3380 bb = b2 * slat * slat 

3390 cc = aa + bb 

3400 dd ~ sqr (cc) 

3410 r = sqr (elevkm * (elevkm + 2. * dd)+(a2 * aa + b2 * bb) / cc) 

3420 cd = (elevkm + dd) / r 

3430 sd = (a2 - b2) / dd * slat * clat / r 

3440 aa = slat 

3450 slat = slat * cd - clat * sd 

3460 clat = clat * cd + aa * sd 

3470 end if 

3480 ratio = erad / r 

3490 aa = sqr (3.) 

3500 p(l) = 2. * slat 

3510 p(2) = 2. * clat 

3520 p(3) = 4.5 * slat * slat - 1.5 

3530 p(4) = 3. * aa * clat * slat 

3540 q(l) = -clat 

3550 q(2) = slat 

3560 q(3) = -3. * clat * slat 

3570 q(4) = aa * (slat * slat - clat * clat) 

3580 for k = 1 to npq ! this loop ends at statement 4110 

3590 if n < m then 

3600 m = 0 

3610 n = n + 1 

3620 rr = ratio A (n + 2) 

3630 fn = n 

3640 end if 

3650 fm = m 

3660 if k >= 5 then 

3670 if m = n then 

3680 aa = sqr (1. - .5 / fm) 

3690 j = k - n - 1 
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3700 p(k) = (1. + 1. / fm) * aa * clat * p(j) 

3710 q(k) = aa * (clat * q(j) + slat / fm * p(j)) 

3720 sl(m) = sl(m-l) * cl(l) + cl(m-l) * sl(l) 

3730 cl(m) = cl(m-l) * cl(l) - sl(m-l) * sl(l) 

3740 else 

3750 aa = sqr (fn * fn - fm * fm) 

3760 bb = sqr ((fn - l.) A 2 - fm * fm) / aa 

3770 cc = (2. * f n - 1.) / aa 

3780 i = k - n 

3790 j = k - 2 * n + 1 

3800 p(k) = (fn + 1.) * (cc * slat / f n * p(i) - bb/ (fn - 1.) * p(j)) 

3810 q(k) = cc* (slat * q(i) - clat / f n * p(i))- bb * q(j) 

3820 end if 

3830 end if 

3840 aa = rr * gh(l) 

3850 daa= rr * dgh(l) 

3860 if m = 0 then 

3870 x = x + aa * q(k) 

3880 dx= dx + daa * q(k) 

3890 z = z - aa * p(k) 

3900 dz = dz - daa*p(k) 

3910 1 = 1 + 1 

3920 else 

3930 bb = rr * gh(l+l) 

3940 dbb = rr * dgh(l+l) 

3950 cc = aa * cl(m) + bb * sl(m) 

3960 dec = daa * cl(m) + dbb * sl(m) 

3970 x = x + cc * q(k) 

3980 dx = dx + dec * q(k) 

3990 z = z - cc * p(k) 

4000 dz = dz - dec * p(k) 

4010 if clat > 0. then 

4020 y = y + (aa * sl(m) - bb * cl(m)) * fm * p(k)/ ((fn + 1.) * clat) 

4030 dy = dy + (daa * sl(m) - dbb * cl(m)) * fm * p(k)/ ((fn + 1.) * clat) 

4040 else 

4050 y = y + (aa * sl(m) - bb * cl(m)) * q(k)* slat 

4060 dy = dy + (daa * sl(m) - dbb * cl(m)) * q(k)* slat 

4070 end if 

4080 1 = 1 + 2 

4090 end if 

4100 m = m + 1 

4110 next k 

4120 x = x + factor * dx 

4130 y = y + factor * dy 

4140 z = z + factor * dz 

4150 aa = x 

4160 x = x*cd + z*sd 

4170 z = z * cd - aa * sd 

4180 h=sqr(x*x+y*y) 

4190 sn =0.0001 

4200 rad=57.29577951 

4210 t=sqr(x*x+y*y+z*z) 

4220 if x=0 and y=0 then dec=0 
4230 if x=0 and y<0 then dec=-90 
4240 if x=0 and y>0 then dec=90 
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4250 

ayx=(y/x) 



4260 

if x>0 and y<0 

then 

dec-(atn(ayx)) 

4270 

if x<0 and y<0 

then 

dec=(atn(ayx)-180) 

4280 

if x<0 and y>0 

then 

dec=(180-atn(ayx)) 

4290 

if x>0 and y>0 

then 

dec=(atn(ayx)) 

4300 

if h=0 and z-0 

then 

dip=0 

4310 

if h=0 and z>0 

then 

dip=-90 

4320 

if h=0 and z<0 

then 

dip=90 

4330 

azh=(z/h) 


4340 

if z=0 then dip 

-0 


4350 

if z>0 then dip 

= abs(atn(azh)) 

4360 

if z<0 then dip 

= -abs(atn(azh)) 

4370 

end sub 




The following is a listing of the contents of file “IGRF95” called for in the subroutine GEOMAGJ. Only 
coefficients through nmax=5 are illustrated here. The sequence is as shown below. 


m 

n 

g 

h 

dg 

dh 

0 

i 

-29682 

0 

17.6 

0 

1 

i 

- 1769 

5310 

13 

-18.3 

0 

2 

- 2197 

0 

-13.2 

0 

1 

2 

3074 

-2356 

3.7 

-15 

2 

2 

1685 

-425 

-.6 

-6.0 

0 

3 

1320 

0 

1.5 

0 

1 

3 

-2268 

-261 

-6.4 

4.1 

2 

3 

1246 

302 

-.2 

2.2 

3 

3 

766 

-406 

-8.1 

-12.1 

0 

4 

941 

0 

.8 

0 

1 

4 

782 

262 

.9 

1.8 

2 

4 

291 

-232 

-6.9 

1.2 

3 

4 

-421 

99 

.5 

2.7 

4 

4 

116 

-301 

-4.6 

-1 

0 

5 

-210 

0 

.8 

0 

1 

5 

352 

44 

.1 

.2 

2 

5 

237 

157 

-1.5 

1.2 

3 

5 

-122 

-152 

-2 

.3 

4 

5 

-167 

-64 

-.1 

1.8 

5 

etc. 

5 

-26 

99 

2.3 

.9 
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